
#ifndef FUNCTIONS_H_
#define FUNCTIONS_H_

#include "VarArray.hpp"
#include "defs.hpp"
#include <Eigen/Dense>
#include <iostream>
#include <vector>

using namespace std;
using namespace Eigen;

void ReadMeshASSCII(const std::string &meshFname, VarArray<Point, 2, 0> &grid);

void InitializeFlow(VarArray<Variable, 2, 0> &flow);

void ThomasAlgorithm(VarArray<double, 2, 0> &coefficient_matrix, double *value,
                     int begin, int end);

void UpdateU(const vector<double> &deltaU, VarArray<Variable, 2, 0> &flow,
             int column);

double CalcNormal(const vector<double> &deltU);

double UpdateVfirstOrder(VarArray<Variable, 2, 0> &flow, const int column,
                         const VarArray<Point, 2, 0> &grid);

void InitializeNextLine(VarArray<Variable, 2, 0> &flow, int column);

double NewtonIterForU(VarArray<Variable, 2, 0> &flow, const int column,
                      const VarArray<Point, 2, 0> &grid, const double re_num);

double UpdateVsecondOrder(VarArray<Variable, 2, 0> &flow, const int column,
                          const VarArray<Point, 2, 0> &grid);

double UpdateVthirdOrder(VarArray<Variable, 2, 0> &flow, const int column,
                         const VarArray<Point, 2, 0> &grid);

double UpdateVthirdOrder2(VarArray<Variable, 2, 0> &flow, const int column,
                          const VarArray<Point, 2, 0> &grid);

class CalcV {
private:
  VarArray<Variable, 2, 0> &flow_;
  const VarArray<Point, 2, 0> &grid_;
  //   int column_{0};
  MatrixXd coef_;
  VectorXd fxk_;
  VectorXd old_v_;
  int size;

public:
  CalcV(VarArray<Variable, 2, 0> &flow, const VarArray<Point, 2, 0> &grid)
      : flow_{flow}, grid_{grid} {
    size = grid.Dim(1) - 2;
    coef_ = MatrixXd::Zero(size, size);
    VectorXd vv(grid.Dim(1) - 2), vvv(grid.Dim(1) - 2);
    fxk_ = vv;
    old_v_ = vvv;
  }
  void CalcCoefficient() {

    coef_(0, 0) = -1;
    coef_(0, 1) = 8;
    coef_(0, 2) = -1;

    coef_(1, 0) = -8;
    coef_(1, 1) = 0;
    coef_(1, 2) = 8;
    coef_(1, 3) = -1;

    for (int i = 2; i < size - 2; ++i) {
      coef_(i, i - 2) = 1;
      coef_(i, i - 1) = -8;
      coef_(i, i) = 0;
      coef_(i, i + 1) = 8;
      coef_(i, i + 2) = -1;
    }
    coef_(size - 2, size - 4) = 1;
    coef_(size - 2, size - 3) = -8;
    coef_(size - 2, size - 2) = 0;
    coef_(size - 2, size - 1) = 7;

    coef_(size - 1, size - 3) = 1;
    coef_(size - 1, size - 2) = -8;
    coef_(size - 1, size - 1) = 7;

    // cout << coef_ << endl;
  }

  void CalcFxk(int column_) {

    double deltx = grid_(column_, 0).x - grid_(column_ - 1, 0).x;
    double rdeltx = 1.0 / deltx;

    const double Jacob1 =
        0.5 * (grid_(column_, 1 + 1).y - grid_(column_, 1 - 1).y);
    fxk_(0) = 24.0 * Jacob1 * rdeltx *
                  (flow_(column_, 1).u - flow_(column_ - 1, 1).u) -
              flow_(column_ - 1, 1).v + 8.0 * flow_(column_ - 1, 2).v -
              flow_(column_ - 1, 3).v;

    const double Jacob2 =
        0.5 * (grid_(column_, 2 + 1).y - grid_(column_, 2 - 1).y);
    fxk_(1) = 12.0 * Jacob2 * rdeltx *
                  (flow_(column_, 2).u - flow_(column_ - 1, 2).u) -
              8.0 * flow_(column_ - 1, 1).v + 8.0 * flow_(column_ - 1, 3).v -
              flow_(column_ - 1, 4).v;

    for (int j = 3; j < size - 1; ++j) {
      double Jacob = 0.5 * (grid_(column_, j + 1).y - grid_(column_, j - 1).y);
      fxk_(j - 1) =
          12.0 * Jacob * rdeltx *
              (flow_(column_, j).u - flow_(column_ - 1, j).u) +
          flow_(column_ - 1, j - 2).v - 8.0 * flow_(column_ - 1, j - 1).v +
          8.0 * flow_(column_ - 1, j + 1).v - flow_(column_ - 1, j + 2).v;
    }

    const double Jacob3 =
        0.5 * (grid_(column_, size - 1 + 1).y - grid_(column_, size - 1 - 1).y);
    fxk_(size - 2) =
        12.0 * Jacob3 * rdeltx *
            (flow_(column_, size - 1).u - flow_(column_ - 1, size - 1).u) +
        flow_(column_ - 1, size - 1 - 2).v -
        8.0 * flow_(column_ - 1, size - 1 - 1).v +
        8.0 * flow_(column_ - 1, size - 1 + 1).v -
        flow_(column_ - 1, size - 1 + 2).v;

    const double Jacob4 =
        0.5 * (grid_(column_, size + 1).y - grid_(column_, size - 1).y);
    fxk_(size - 1) = 12.0 * Jacob4 * rdeltx *
                         (flow_(column_, size).u - flow_(column_ - 1, size).u) +
                     flow_(column_ - 1, size - 2).v -
                     8.0 * flow_(column_ - 1, size - 1).v +
                     8.0 * flow_(column_ - 1, size + 1).v;

    fxk_ *= -1.0;
  }

  double UpdateVFourthOrder(int column_) {

    int npoint_in_ydir = grid_.Dim(1);

    // vector<double> old_v(npoint_in_ydir - 2);
    for (int j = 1; j < npoint_in_ydir - 1; ++j)
      old_v_(j - 1) = flow_(column_, j).v;

    // Vector3f x = A.colPivHouseholderQr().solve(b);
    VectorXd new_v = coef_.colPivHouseholderQr().solve(fxk_);

    double res_v{0.0};
    for (int j = 0; j < npoint_in_ydir - 2; ++j) {
      res_v += (new_v(j) - old_v_(j)) * (new_v(j) - old_v_(j));
      flow_(column_, j + 1).v = new_v(j);
    }

    // for (int j = 0; j < npoint_in_ydir; ++j)
    //   cout << flow_(column_, j).v << endl;
    // cout << "end v" << endl;
    // exit(EXIT_FAILURE);
    return sqrt(res_v);
  }
};

#endif